function main_continuation_dro_
clc;
close all;

% mu = 3.04036e-6; % parameter of system (Sun-Earth)
mu = 1.21506683e-2; % parameter of system (Earth-Moon)
RelTol = 1e-13; % tolerances of integration
AbsTol = 1e-13; % tolerances of integration
MaxNum = 25; % maximum number of iterations
Maxtf = 10; % maximum time of flight

auxdata.mu = mu;
auxdata.RelTol = RelTol;
auxdata.AbsTol = AbsTol;
auxdata.MaxNum = MaxNum;
auxdata.Maxtf = Maxtf;

options = odeset('RelTol',RelTol,'AbsTol',AbsTol,'Events',@event_);

% 延拓参数
    continuation_num = 0;
    step_size = 0.01; % 延拓步长
    max_continuation = 100; % 最大延拓次数

% dro初值
x0 = [-7e4/384400+(1-mu) ; 0 ; 0 ; 0 ; 0.5 ; 0]; %月心

% 微分修正后的dro精确初值
[x0_new, period] = dro_(x0 , auxdata);

figure('position' , [100 , 100 , 400 , 280]) % 绘制区域的位置和大小
x0_new_continuation = x0_new;
% x作为延拓参数
while (continuation_num <= max_continuation)
      % 更新步长
      if continuation_num > 20
         step_size = 0.0025; % 减小步长
      elseif continuation_num > 25
          step_size = 0.001;
      elseif continuation_num > 48
          step_size = 0.00005;
      end

      % 更新初始猜测值
      x0_new_ = x0_new_continuation;
      x0_new_(1) = x0_new_(1) - step_size
      
      % 使用微分修正得到新的 DRO 初值
      [x0_new_continuation, period] = continuation_dro_(x0_new_, auxdata);

      %利用修正后的延拓初值求解crtbp
      [~ , xx0_new_continuation] = ode45(@fun_crtbp_ , [0 , 100] , x0_new_continuation , options , auxdata); 

      % 画dro轨道
      plot3(xx0_new_continuation(: , 1) , xx0_new_continuation(: , 2) , xx0_new_continuation(: , 3) , 'k' , 'LineWidth' , 2); hold on;
      xlabel('x/LU'); ylabel('y/LU'); zlabel('z/LU');
      axis equal;
      axis([-1 2.2 , -1.8 1.8 -0.1 0.1])
      hold on

      continuation_num = continuation_num + 1

end

% 积分终止条件，即积分到x-z平面终止
    function [value,isterminal,direction] = event_(t, xx, auxdata)
    value = xx(2); % 检测 x(2) = 0 
    isterminal = 1; % 终止积分
    direction = 1; % 正向穿过
    end

end
